#ifndef GUARD_Complex_Polynomial
#define GUARD_Complex_Polynomial

// Complex_Polynomial.h header file
// Licensed under GPLv2
// Terms available at http://www.gnu.org/licenses/gpl-2.0.txt

#include <vector>
#include <complex>

// Complex_Polynomial stores monic polynomials based on their roots,
// for ease of computation. This does limit the ability of the code to
// generate fractals with more interesting attributes, but in this
// context I felt that that limitation was worthwhile. After all, this
// isn't being done in a CAS where factoring is an option (unlike most
// of the Mathematica code I've seen for this.)

struct Complex_Polynomial {
  std::vector<std::complex<double> > roots;

  std::complex<double> evalPt( std::complex<double> pt ) {
    //Gotta make sure I evaluate the whole polynomial
    std::vector<std::complex<double> >::size_type i = 0;
    //Sensible starting solution
    std::complex<double> val ( 1 , 1 );

    for ( i = 0 ; i < roots.size() ; i++ ) {
      //Multiplying each factor by the result so far
      val *= ( pt - roots[i] );
    }
    return val;
  }

  std::complex<double> derivEval( std::complex<double> pt ) {
    //I'll need two indices to track which I'm evaluating and which
    //I'm avoiding.
    std::vector<std::complex<double> >::size_type i,j,k;
    if ( roots.size() == 0 ) {
      return std::complex<double> ( 0 , 0 );
    }

    //This is the first stage of derivative evaluation.
    std::vector<std::complex<double> > stages (roots.size());

    // Since the derivative of a polynomial by its roots is defined by
    // the sum of the polynomials formed by each set of roots less
    // one, this is the fastest route with my arrangement.
    for ( k = 0 ; k < stages.size() ; k++ ) {
      // Filling all the roots with initial solutions for multiplication.
      stages[k] = std::complex<double> ( 1 , 1 ) ;
    }

    // Here's the implementation, and it relies on skipping certain
    // elements. This works by evaluating every factor but the ith,
    // and inputting that into the ith position of the stages vector.

    for ( i = 0 ; i < roots.size() ; i ++ ) {
      for ( j = 0 ; j < roots.size() ; j++ ) {
        if ( j != i ) {
          stages[i] *= ( pt - roots[j] ) ;
        }
      }
    }
    // Now that all the evaluations have occurred, the summation can happen!
    std::complex<double> result ( 0 , 0 );
    for ( k = 0 ; k < stages.size() ; k++ ) {
      //Sums all the components of the derivative at the value.
      result += stages[k];
    }
    return result;
  }
};

std::complex<double> evalPt( std::complex<double> pt );
std::complex<double> derivEval( std::complex<double> pt );

#endif
